#include <test_csm.hpp>

#include <life/lifediscr/mesh.hpp>






namespace Life
{


template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void
Test_Schwarz<Dim, Order, Entity>::run()
{   
    using namespace Life::vf;

    //----------------------------------------------------------------------------------------------------------//
    
    if ( this->vm().count( "help" ) )
        {
            std::cout << this->optionsDescription() << "\n";
            return;
        }

    //----------------------------------------------------------------------------------------------------------//
    
    mesh_ptrtype mesh = createMesh( meshSize );
    mesh_ptrtype mesh1 = createSubMesh(mesh,mesh->markerName("Physique1"));
    mesh_ptrtype mesh2 = createSubMesh(mesh,mesh->markerName("Physique2"));

    space_ptrtype Xh1 = space_type::New( mesh1 );
    space_ptrtype Xh2 = space_type::New( mesh2 );

    element_type U1( Xh1, "u1" );
    element_type U2( Xh2, "u2" );

    AUTO( SolAnalytique, cos(M_PI*Px()/2.0)*cos(M_PI*Py()/2.0) );
    U1 = project( Xh1, elements(Xh1->mesh()), SolAnalytique );

    //double __integrate = integrate(elements(mesh1),
    //idv(U1) ).evaluate()(0,0);
    double __integrate_bd_m1 = integrate(boundaryfaces(mesh1),
				   idv(U1) ).evaluate()(0,0);

    double __integrate_m = integrate(markedfaces(mesh,mesh->markerName("Interface")),
    				   idv(U1) ).evaluate()(0,0);
    double __integrate_m1 = integrate(markedfaces(mesh1,mesh->markerName("Interface")),
    				   idv(U1) ).evaluate()(0,0);
    std::cout<<"\nIntergate bd m1=" << __integrate_bd_m1<<"\n";
    std::cout<<"\nIntergate m=" << __integrate_m <<" et Integrate m1=" << __integrate_m1 <<"\n";

    
    std::cout << "\nMESH1 boundary";
    typename MeshTraits<mesh_type>::location_face_const_iterator itf2, enf2;
    boost::tie( boost::tuples::ignore, itf2, enf2 ) = boundaryfaces( mesh1 );
    for ( ; itf2!=enf2 ; ++itf2)
      {
	//std:: cout << "\nface p0 id : " <<itf2->point(0).showMe();//id();
	//std:: cout << "\nface p1 id : " <<itf2->point(1).showMe();//id();
	std:: cout << "\nface p0 id : " <<itf2->point(0).id();
	std:: cout << "\nface p1 id : " <<itf2->point(1).id();

      }



    std::cout << "\nMESH";
    typename MeshTraits<mesh_type>::marker_face_const_iterator itf, enf;
    boost::tie( boost::tuples::ignore, itf, enf ) = markedfaces( mesh, mesh->markerName("Interface"));
    for ( ; itf!=enf ; ++itf)
      {
	std:: cout << "\nface p0 id : " <<itf->point(0).showMe();//id();
	std:: cout << "\nface p1 id : " <<itf->point(1).showMe();//id();
      }
    std::cout << "\nMESH1";
    //typename MeshTraits<mesh_type>::marker_face_const_iterator itf, enf;
    boost::tie( boost::tuples::ignore, itf, enf ) = markedfaces( mesh1, mesh->markerName("Interface"));
    for ( ; itf!=enf ; ++itf)
      {
	std:: cout << "\nface p0 id : " <<itf->point(0).showMe();//id();
	std:: cout << "\nface p1 id : " <<itf->point(1).showMe();//id();
      }


    //    test_data_maillage();

    // test_interpolation();

    // test_composantes_champs();
    

    typename MeshTraits<mesh_type>::marker_element_const_iterator it, en;
    boost::tie( boost::tuples::ignore, it, en ) = markedelements( mesh, mesh->markerName("Physique1"));

    //      for ( ; it!=en;++it )
      {
	std::cout<<"\n----------Element id " << it->id() << "----------";
	std::cout<<"\nPOINTS:\n";
	for (unsigned int n=0; n < it->nPoints(); n++) 
	  {
	    std::cout<<it->point(n).id()<<"\n";//.showMe();
	  }

	//std::cout<<"\nFACES:"<<it->nFaces()<<"\n";
	std::cout<<"\nFACES:"<<it->nNeighbors()<<"\n";std::cout<<"\nfaceid:"<<it->face(0).id()<<"\n";
	for (unsigned int s=0; s<it->nNeighbors(); s++)
	//for (unsigned int s=0; s< 3/*it->nFaces()*/; s++)
	  {                            
    
	    std::cout<<"\nface num0="<<boost::get<0>( mesh->localFaceId( it->id(), s ) );
	    //std::cout<<"\n  face num1="<<boost::get<1>( mesh->localFaceId( it->id(), s ) )<<"\n";

	    typename mesh_type::face_type new_face = mesh->face( boost::get<0>( mesh->localFaceId( it->id(), s ) ) );
	    //typename mesh_type::face_type new_face = mesh->face( 5 );
	    for ( uint16_type p = 0;p < 2/*new_face.nPoints()*/; ++p )
	      {
		//std::cout<<it->face(s).point(p).id()<<"\n";//.showMe();
		//std::cout<<new_face.point(p).id()<<"\n";//showMe();
		//mesh->face( boost::get<1>( mesh->localFaceId( it->id(), s ) )).point(p).showMe();
		std::cout<<mesh->face(11).point(p).id()<<"\n";
	      }
	  }
	       }





} // Test_Schwarz::run


    //----------------------------------------------------------------------------------------------------------//
    //----------------------------------------------------------------------------------------------------------//
    //----------------------------------------------------------------------------------------------------------//
    //----------------------------------------------------------------------------------------------------------//
    //----------------------------------------------------------------------------------------------------------//


template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void 
Test_Schwarz<Dim, Order, Entity>::test_data_maillage() {

    // First we create the mesh
    mesh_ptrtype mesh = createMesh( meshSize );
    mesh_ptrtype mesh1 = createSubMesh(mesh,mesh->markerName("Physique1"));
    mesh_ptrtype mesh2 = createSubMesh(mesh,mesh->markerName("Physique2"));
    
    //----------------------------------------------------------------------------------------------------------//

    typedef typename mesh_type::element_const_iterator element_const_iterator;
    typedef typename mesh_type::element_type mesh_element_type;

    std::cout<<"\nElement marker :\nMesh\n";
    for ( element_const_iterator elt = mesh->beginElement();  elt != mesh->endElement(); ++elt ) {
    std::cout<<" "<<elt->marker().value(); }

    std::cout<<"\nMesh1\n";
    for ( element_const_iterator elt = mesh1->beginElement();  elt != mesh1->endElement(); ++elt ) {
    std::cout<<" "<<elt->marker().value(); }

    std::cout<<"\nMesh2\n";
    for ( element_const_iterator elt = mesh2->beginElement();  elt != mesh2->endElement(); ++elt ) {
    std::cout<<" "<<elt->marker().value(); }

    //----------------------------------------------------------------------------------------------------------//

    typedef typename mesh_type::face_const_iterator face_const_iterator;

    std::cout<<"\n\nFace marker :\nMesh\n";
    for ( face_const_iterator elt = mesh->beginFace();  elt != mesh->endFace(); ++elt ) {
    std::cout<<" "<<elt->marker().value(); }

    std::cout<<"\nMesh1\n";
    for ( face_const_iterator elt = mesh1->beginFace();  elt != mesh1->endFace(); ++elt ) {
    std::cout<<" "<<elt->marker().value(); }

    std::cout<<"\nMesh2\n";
    for ( face_const_iterator elt = mesh2->beginFace();  elt != mesh2->endFace(); ++elt ) {
    std::cout<<" "<<elt->marker().value(); }

    //----------------------------------------------------------------------------------------------------------//

    typedef typename mesh_type::point_const_iterator point_const_iterator;

    std::cout<<"\n\nPoint marker :\nMesh\n";
    for ( point_const_iterator elt = mesh->beginPoint();  elt != mesh->endPoint(); ++elt ) {
    std::cout<<" "<<elt->marker().value(); }

    std::cout<<"\nMesh1\n";
    for ( point_const_iterator elt = mesh1->beginPoint();  elt != mesh1->endPoint(); ++elt ) {
    std::cout<<" "<<elt->marker().value(); }

    std::cout<<"\nMesh2\n";
    for ( point_const_iterator elt = mesh2->beginPoint();  elt != mesh2->endPoint(); ++elt ) {
    std::cout<<" "<<elt->marker().value(); }

    //----------------------------------------------------------------------------------------------------------//

    typedef typename mesh_type::super_faces::location_face_iterator location_face_iterator;

    std::cout<<"\n\nBoundary Face marker :\nMesh\n";
    for ( location_face_iterator elt = mesh->facesOnBoundary().first;  elt != mesh->facesOnBoundary().second; ++elt ) {
    std::cout<<" "<<elt->marker().value(); }


  //  typedef typename mesh_type::super_faces::location_face_iterator location_face_iterator;

    std::cout<<"\n\nInternal Face marker :\nMesh\n";
    for ( location_face_iterator elt = mesh->internalFaces().first;  elt != mesh->internalFaces().second; ++elt ) {
    std::cout<<" "<<elt->marker().value(); }

    //----------------------------------------------------------------------------------------------------------//

    std::cout<<"\n\n\n";
    for ( element_const_iterator elt = mesh->beginElement();  elt != mesh->endElement(); ++elt ) {
        std::cout<<"\n";
        //for (unsigned int n=0; n < elt->nFaces(); n++) {
        //  for (unsigned int s=0; s<elt->nNeighbors(); s++) {
        // std::cout<< (mesh->face( boost::get<0>( mesh->localFaceId( elt->id(), s ) ) ) ).marker().value() <<" " ; }
            //  (mesh->face( boost::get<0>( mesh->localFaceId( elt->id(), s ) ) ) ) ; 
        //    std::cout<<mesh->face( 0 ).marker().value()<<" " ;
          //  std::cout<<mesh->face( 1 ).marker().value()<<" " ; 
           // boost::get<0>( mesh->localFaceId( elt->id(), s ) );
        std::cout<<" "<<elt->nFaces();
        std::cout<<" "<<elt->nEdges();
        std::cout<<" "<<elt->nVertices();
        std::cout<<" "<<elt->nPoints()<<" /// ";
        //elt->face(n);
        //mesh->localFaceId( elt->id(), s );
        //  }
    }



}


template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void 
Test_Schwarz<Dim, Order, Entity>::test_interpolation() {

using namespace Life::vf;


    mesh_ptrtype mesh1 = createMesh1( meshSize );
    mesh_ptrtype mesh2 = createMesh2( meshSize );

    space_ptrtype Xh1 = space_type::New( mesh1 );
    space_ptrtype Xh2 = space_type::New( mesh2 );

    //----------------------------------------------------------------------------------------------------------//
  
    element_type U1( Xh1, "u1" );element_type U1b( Xh1, "u1b" );element_type U1ex( Xh1, "u1ex" );
    element_type V1( Xh1, "v1" );
    element_type U2( Xh2, "u2" );element_type U2b( Xh2, "u2b" );element_type U2ex( Xh2, "u2ex" );
    element_type V2( Xh2, "v2" );

    //----------------------------------------------------------------------------------------------------------//

    AUTO( SolAnalytique, cos(M_PI*Px()/2.0)*cos(M_PI*Py()/2.0) );
    U1ex = project( Xh1, elements(Xh1->mesh()), SolAnalytique );
    U2ex = project( Xh2, elements(Xh2->mesh()), SolAnalytique );

    U1 = project( Xh1, elements(Xh1->mesh()), SolAnalytique );

    interpolate<space_type,element_type>( Xh2, U1, U2 );
    U2.updateGlobalValues();//???interet de cette fonction???
    
    std::cout<<"\nDoit valoir 0 (idv) :"
             <<math::sqrt( integrate(markedfaces(mesh2,mesh2->markerName("Interface")), _Q<2*Order>(),
                                    (idv(U2)-idv(U2ex))*(idv(U2)-idv(U2ex)) ).evaluate()(0,0) ) <<"\n";
    //    (U2(0)-U2(0))*(idv(U2)-idv(U2ex)) ).evaluate()(0,0) ) <<"\n";

    std::cout<<"\nDoit valoir 0 (dxv):"
             <<math::sqrt( integrate(markedfaces(mesh2,mesh2->markerName("Interface")), _Q<1*Order>(),
    //                              (gradv(U2)-gradv(U2ex))*trans(gradv(U2)-gradv(U2ex)) ).evaluate()(0,0) ) <<"\n";
                                    (dxv(U2)-dxv(U2ex))*trans(dxv(U2)-dxv(U2ex)) ).evaluate()(0,0) ) <<"\n";

    //----------------------------------------------------------------------------------------------------------//

    typedef typename mesh_type::element_type geoelement_type;
    // geometric mapping context
    typedef typename mesh_type::gm_type gm_type;
    typedef boost::shared_ptr<gm_type> gm_ptrtype;
    typedef typename gm_type::template Context<vm::POINT, geoelement_type> gmc_type;
    typedef boost::shared_ptr<gmc_type> gmc_ptrtype;

    typedef typename space_type::basis_type space_basis_type;
    space_basis_type const* __basis = Xh1->basis().get();

    gm_ptrtype __gm = Xh1->gm();
    typedef typename gm_type::precompute_ptrtype geopc_ptrtype;
    typedef typename gm_type::precompute_type geopc_type;
    geopc_ptrtype __geopc( new geopc_type( __gm, __basis->dual().points() ) );


    typedef typename element_type::id_type id_type;

    typedef typename mesh_type::element_const_iterator element_const_iterator;
    typedef typename mesh_type::element_type mesh_element_type;

    for ( element_const_iterator elt = mesh1->beginElement();  elt != mesh1->endElement(); ++elt ) {

        gmc_ptrtype __c( new gmc_type( __gm, *elt, __geopc ) );
        typename element_type::pc_type pc( U1.functionSpace()->fe(), __c->xRefs() );
        typename element_type::id_type interpfunc( U1.id( *__c, pc ) );

        cout<<"\n ooo :"<<interpfunc.M_id.size()<<" "<<interpfunc.M_id.num_dimensions()
            <<" "<<interpfunc.M_id.num_elements()<<" "<<interpfunc(0,0,0)<<" "
            <<interpfunc(0,0,1)<<" "<<interpfunc(0,0,2)<<" "<<interpfunc(0,0,3)<<" " ;

    }
}


template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
void 
Test_Schwarz<Dim, Order, Entity>::test_composantes_champs() {

    using namespace Life::vf;

    mesh_ptrtype mesh = createMesh1( meshSize );

    //Test composante d'un champ vectoriel
    space_ptrtype2 XhV = space_type2::New( mesh );
    element_type2 UVec( XhV, "uvec" );

    //AUTO( FonctionVec, vec( cos(M_PI*Px()/2.0), cos(M_PI*Py()/2.0) ) );
    //    AUTO( FonctionVec, vec( Px(), Py() ) );
    AUTO( FonctionVec, vec( cst(1.0), cst(-1.0) ) );

    typedef typename element_type2::component_type component_type;

    UVec=project( XhV, elements(XhV->mesh()), FonctionVec );
    UVec.updateGlobalValues();

    cout<<"\natat "<< UVec.start()+Y*XhV->nDofPerComponent();

    component_type uXXX = UVec.template comp<X>(boost::mpl::bool_<false>());
    component_type uYYY = UVec.template comp<Y>(boost::mpl::bool_<false>());


    std::cout<<"\nDoit valoir 0  Comp:"
             <<math::sqrt( integrate( elements(mesh), _Q<5*Order>(), 
                           trans( idv(UVec) - vec(idv(uXXX),idv(uYYY)) ) * ( idv(UVec) - vec(idv(uXXX),idv(uYYY)) ) ).evaluate()(0,0) ) <<"\n";

    std::cout<<"\nDoit valoir 0  Comp1:"
             <<math::sqrt( integrate( elements(mesh), _Q<5*Order>(), 
                           trans( trans(idv(UVec))*oneX() - idv(uXXX) ) * ( trans(idv(UVec))*oneX() - idv(uXXX) ) ).evaluate()(0,0) ) <<"\n";

    std::cout<<"\nDoit valoir 0  Comp2:"
             <<math::sqrt( integrate( elements(mesh), _Q<5*Order>(), 
                           trans( /*trans(idv(UVec))*oneY() -*/ idv(uYYY) ) * ( /*trans(idv(UVec))*oneY() -*/ idv(uYYY) ) ).evaluate()(0,0) ) <<"\n";

    std::cout<<"\nDoit valoir 0  Comp2:"
             <<math::sqrt( integrate( elements(mesh), _Q<5*Order>(), 
                             ( trans(idv(UVec))*oneY() ) * ( trans(idv(UVec))*oneY() ) ).evaluate()(0,0) ) <<"\n";


    std::cout<<"\nNbdof de UVec : "<< UVec.nDof();
    std::cout<<"\nNbdof de uXXX : "<< uXXX.nDof();
    std::cout<<"\nNbdof de uYYY : "<< uYYY.nDof();

    for( unsigned int i = 0; i < uXXX.nDof(); ++i ) { std::cout<<"\n"<<uXXX(i)<<" "<<UVec(i);}
    cout<<"\n----------------------------------------------";
    for( unsigned int i = 0; i < uYYY.nDof(); ++i ) { std::cout<<"\n"<<uYYY(i)<<" "<<UVec(uXXX.nDof() + i);}





   //----------------------------------------------------------------------------------------------------------//

    typedef typename mesh_type::element_type geoelement_type;
    // geometric mapping context
    typedef typename mesh_type::gm_type gm_type;
    typedef boost::shared_ptr<gm_type> gm_ptrtype;
    typedef typename gm_type::template Context<vm::POINT|vm::JACOBIAN|vm::HESSIAN|vm::KB, geoelement_type> gmc_type;
//    typedef typename gm_type::template Context<vm::POINT, geoelement_type> gmc_type;
    typedef boost::shared_ptr<gmc_type> gmc_ptrtype;

    typedef typename space_type2::basis_type space_basis_type;
    space_basis_type const* __basis = XhV->basis().get();
     typedef typename space_type2::component_basis_type component_basis_type;
    component_basis_type const* __basis_comp =  UVec.compSpace()->basis().get();

    gm_ptrtype __gm = XhV->gm();
    typedef typename gm_type::precompute_ptrtype geopc_ptrtype;
    typedef typename gm_type::precompute_type geopc_type;
    geopc_ptrtype __geopc( new geopc_type( __gm, __basis->dual().points() ) );
    gm_ptrtype __gm_comp = XhV->gm();//UVec.compSpace()->gm();
    geopc_ptrtype __geopc_comp( new geopc_type( __gm_comp, __basis_comp->dual().points() ) );


    typedef typename element_type2::id_type id_type;
    typedef typename component_type::id_type id_type_comp;

    typedef typename mesh_type::element_const_iterator element_const_iterator;
    typedef typename mesh_type::element_type mesh_element_type;

    cout<<"\n\n\n";
    for ( element_const_iterator elt = mesh->beginElement();  elt != mesh->endElement(); ++elt ) {

        gmc_ptrtype __c( new gmc_type( __gm, *elt, __geopc ) );
        typename element_type2::pc_type pc( UVec.functionSpace()->fe(), __c->xRefs() );
        typename element_type2::id_type interpfunc( UVec.id( *__c, pc ) );

        gmc_ptrtype __c_comp( new gmc_type( __gm/*_comp*/, *elt, __geopc/*_comp*/ ) );
        typename component_type::pc_type pc_comp( UVec.compSpace()->fe(), __c/*_comp*/->xRefs() );
        typename component_type::id_type interpfunc_comp( uXXX.id( *__c_comp, pc_comp ) );

      /*  cout<<"\n ooo :"<<interpfunc.M_id.size()<<" "<<interpfunc.M_id.num_dimensions()
            <<" "<<interpfunc.M_id.num_elements()<<" "<<interpfunc(0,0,0)-interpfunc_comp(0,0,0)<<" "
            <<interpfunc(0,0,1)<<" "<<interpfunc(0,0,2)<<" "<<interpfunc(0,0,3)<<" " ;*/
        for (int i=0;i<6;++i) {        
        cout<<"\n"<<interpfunc(0,0,i)<<" "<<interpfunc_comp(0,0,i);}
        cout<<"\n----------------------------------------------";
    }

    //----------------------------------------------------------------------------------------------------------//

    for( unsigned int i = 0; i < uXXX.nDof(); ++i ) { std::cout<<"\n"<<uXXX.globalValue(i)<<" "<<UVec.globalValue(i);}
        cout<<"\n----------------------------------------------";
    //    for( unsigned int i = 0; i < uYYY.nDof(); ++i ) { std::cout<<"\n"<<uYYY(i)<<" "<<UVec(uXXX.nDof() + i);}
    cout<<"\n"<<uXXX.nComponents1;
    cout<<"\n----------------------------------------------\n";

    //----------------------------------------------------------------------------------------------------------//

/*    for ( element_const_iterator elt = mesh->beginElement();  elt != mesh->endElement(); ++elt ) {
        //gmc_ptrtype context( new gmc_type( __gm_comp, *elt, __geopc_comp ) );
        gmc_ptrtype context( new gmc_type( __gm, *elt, __geopc ) );
        for( int l = 0; l < space_basis_type::nDof; ++l ) {

            for( typename element_type2::array_type::index c1 = 0; c1 < UVec.nComponents1; ++c1 ) {

             typename element_type2::array_type::index ldof = space_basis_type::nDof*c1+l;
             size_type gdof = boost::get<0>( UVec.functionSpace()->dof()->localToGlobal( context->id(), l, c1 ) );
                   
         //    double v_ = this->globalValue( gdof );
                cout<<"\n"<<UVec.globalValue( gdof );


        
      }
    }
  }
*/
    //----------------------------------------------------------------------------------------------------------//

    for ( element_const_iterator elt = mesh->beginElement();  elt != mesh->endElement(); ++elt ) {
        gmc_ptrtype context( new gmc_type( __gm_comp, *elt, __geopc_comp ) );
        for( int l = 0; l < component_basis_type::nDof; ++l ) {

            for( typename component_type::array_type::index c1 = 0; c1 < uXXX.nComponents1; ++c1 ) {

                //typename component_type::array_type::index ldof = component_basis_type;::nDof*c1+l;
//                size_type gdof = boost::get<0>( uXXX.functionSpace()->dof()->localToGlobal( context->id(), l, c1 ) );
                size_type gdof = boost::get<0>( uXXX.functionSpace()->dof()->localToGlobal( elt->id(), l, c1 ) );
                   
                //double v_ = this->globalValue( gdof );
                cout<<"\n"<<uXXX.globalValue( gdof )<<"   "<<uXXX( gdof );
            }
        }
    }



}



} // Life


// instantiation
template class Life::Test_Schwarz<2,2,Life::Simplex>;
//template class Life::Test_Schwarz<3,2,Life::Simplex>;
